options(future.globals.maxSize= 891289600)
### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
install.packages("Seurat")
if (!requireNamespace("dplyr", quietly = TRUE))
install.packages("dplyr")
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
devtools::install_github("jinworks/CellChat")
if (!requireNamespace("ComplexHeatmap", quietly = TRUE))
BiocManager::install("ComplexHeatmap")
### Load all the necessary libraries
library(Seurat)
library(dplyr)
library(CellChat)
library(ComplexHeatmap)6 - Cell-Cell Communication analysis
Introduction
In this notebooks we are going to carry out a cell-cell communication analysis using CellChat. You can see the GitHub repository here, the CellChat vignette here and the differential communication analysis vignette here. The goal of the notebook is to identify which communication pathways are altered between flu and covid infection!
Some associated literature which is a must read are:
Key Takeaways
The Ligand-Receptor (L-R) database used gather the collective prior knowledge and has a great impact on the results obtained.
Different CCC tools have varying assumptions, therefore, the tool of choise will also have a major impact on the results.
CellChatandCellPhoneDBmake a point of modelling CCC events taking into account heteromeric complexes. This ensures all the subunits of a protein complex are expressed to consider a cell-cell interaction feasible. This assumption reduces false positive predictions.CellChat, additionally, accounts for interaction mediator proteins such as agonists.
Broadly, CCC tools are generally able to capture relevant biological signals. However, predicted interactions tend to have false positives, if available leveraging information from additional modalities and analyses could help to refine the predictions.
CCC inference from scRNAseq data makes the assumption that gene expression is a proxy for protein levels. Moreover, they don’t (and can’t) account for other intermediate steps between translation and protein function such as post-translational modifications, secretion, diffusion…
Library
Load data
We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.
# Download the data in data/ directory
# download.file(
# url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
# destfile = "../data/workshop-data.rds",
# method = "wget",
# extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds
se <- readRDS("../data/workshop-data.rds")Analysis
Convert ENSEMBL IDs to Gene Symbols
Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:
head(rownames(se))[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
[5] "ENSG00000000460" "ENSG00000000938"Convert to gene symbols
gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")
symbol_id <- data.frame(index = rownames(se)) %>%
left_join(gene_df) %>%
pull(feature_name)
# re-create seurat object
mtx <- se@assays$RNA@data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)# symbol_id <- mapIds(
# org.Hs.eg.db,
# keys = rownames(se),
# column = "SYMBOL",
# keytype = "ENSEMBL",
# multiVals = "first")
#
# # df <- data.frame(symbol = symbol_id, ensembl = names(symbol_id))
# all(rownames(se) == names(symbol_id))
#
# # re-create seurat object
# mtx <- se@assays$RNA@data
# rownames(mtx) <- symbol_id
#
# # Remove NAs
# sum(is.na(rownames(mtx)))
# dim(mtx)
# mtx <- mtx[!is.na(rownames(mtx)), ]
# dim(mtx)
#
# rownames(mtx) <- make.unique(rownames(mtx))
# se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
#
# sum(is.na(rownames(se)))Split dataset
# First we will remove the Uncategorized1/2 populations as we have seen that they are doublets
se <- se[, ! se$Celltype %in% c("Uncategorized1", "Uncategorized2")]
# Create samples so the CellChat knows which cells come from each sample
se$samples <- se$`Sample ID`
# Convert to character to drop inexisting factors
se$Celltype <- factor(as.character(se$Celltype))# Subset to cells coming form COVID-19 patients
covid <- se[, se$disease == "COVID-19"]
# Subset to cells coming form flu patients
flu <- se[, se$disease == "influenza"]Let’s visualize the cell type populations between these 2 datasets:
bind_rows(
data.frame(prop.table(table(covid$Celltype)), disease = "covid"),
data.frame(prop.table(table(flu$Celltype)), disease = "flu")) %>%
ggplot(aes(x = Var1, y = disease, fill = Freq, label = round(Freq, 2))) +
geom_tile(color = "lightgrey") +
geom_text(color = "lightgrey") +
scale_fill_viridis_c(option = "magma") +
labs(title = "Cell type proportions", x = "", y = "") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
Carry out pre-processing steps on covid data and create the cellchat object
# Normalize data
covid <- NormalizeData(covid, verbose = FALSE)
cc_covid <- createCellChat(
object = covid,
group.by = "Celltype",
assay = "RNA",
do.sparse = TRUE)[1] "Create a CellChat object from a Seurat object"
The `meta.data` slot in the Seurat object is used as cell meta information
Set cell identities for the new CellChat object
The cell groups used for CellChat analysis are B cell, IgG-, B cell, IgG+, CD4, EM-like, CD4, non-EM-like, CD8, EM-like, CD8, non-EM-like, classical Monocyte, DC, intermediate Monocyte, NK cell, nonclassical Monocyte, Platelet, RBC Same for the flu
# Normalize data
flu <- NormalizeData(flu, verbose = FALSE)
cc_flu <- createCellChat(
object = flu,
group.by = "Celltype",
assay = "RNA",
do.sparse = TRUE)[1] "Create a CellChat object from a Seurat object"
The `meta.data` slot in the Seurat object is used as cell meta information
Set cell identities for the new CellChat object
The cell groups used for CellChat analysis are B cell, IgG-, B cell, IgG+, CD4, EM-like, CD4, non-EM-like, CD8, EM-like, CD8, non-EM-like, classical Monocyte, DC, intermediate Monocyte, NK cell, nonclassical Monocyte, Platelet, RBC CellChat CCC analysis
Cellchat intuition
How is the communication probability computed? See the methods section in the CellChat paper.
The communication probability Pi,j from cell groups i to j for a particular ligand-receptor pair k was modeled by
\[ P_{i,j}^k = \frac{{L_iR_j}}{{K_h + L_iR_j}} \times \left( {1 + \frac{{AG_i}}{{K_h + AG_i}}} \right) \cdot \left( {1 + \frac{{AG_j}}{{K_h + AG_j}}} \right) \\ \times \frac{{K_h}}{{K_h + AN_i}} \cdot \frac{{K_h}}{{K_h + AN_j}} \times \frac{{n_in_j}}{{n^2}} \]
\[ {L_i = \root {{m1}} \of {{L_{i,1} \cdots L_{i,m1}}},\,R_j = \root {{m2}} \of {{R_{j,1} \cdots R_{j,m2}}} \cdot \frac{{1 + RA_j}}{{1 + RI_j}}} \]
Li and Rj represent the expression level of ligand L and receptor R in cell group i and cell group j, respectively
Kh whose default value was set to be 0.5
AG is the average expression of the L-R agonists
CellChat Database
Define the use of the human L-R database. The database used is a key decision in the analysis workflow since it will define the universe of Ligand-Receptor interactions we are going to be looking for. We want it to be as comprehensive as possible while simultaneously only including curated communication partners.
We personally like the CellChatDB database because:
Groups LR interactions into functional pathways.
Takes into account multimeric complexes as well as other cofactors, such as agonists and antagonists present in the dataset.
LR interactions come from KEGG Pathway database and manually curated from primary references.
We can take a glimpse at all of this information by looking at the database loaded into R!
CellChatDB <- CellChatDB.human
showDatabaseCategory(CellChatDB)
# glimpse(CellChatDB)CellChatDB$interaction %>%
dplyr::select(interaction_name, pathway_name, ligand, receptor, agonist, evidence, annotation) %>%
head() %>%
as_tibble()# A tibble: 6 × 7
interaction_name pathway_name ligand receptor agonist evidence annotation
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 TGFB1_TGFBR1_TGFBR2 TGFb TGFB1 TGFbR1_R2 TGFb a… KEGG: h… Secreted …
2 TGFB2_TGFBR1_TGFBR2 TGFb TGFB2 TGFbR1_R2 TGFb a… KEGG: h… Secreted …
3 TGFB3_TGFBR1_TGFBR2 TGFb TGFB3 TGFbR1_R2 TGFb a… KEGG: h… Secreted …
4 TGFB1_ACVR1B_TGFBR2 TGFb TGFB1 ACVR1B_TG… TGFb a… PMID: 2… Secreted …
5 TGFB1_ACVR1C_TGFBR2 TGFb TGFB1 ACVR1C_TG… TGFb a… PMID: 2… Secreted …
6 TGFB2_ACVR1B_TGFBR2 TGFb TGFB2 ACVR1B_TG… TGFb a… PMID: 2… Secreted …CellChatDB$interaction %>%
dplyr::select(interaction_name, pathway_name, ligand, receptor, agonist, evidence, annotation) %>%
tail() %>%
as_tibble()# A tibble: 6 × 7
interaction_name pathway_name ligand receptor agonist evidence annotation
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Testosterone-Testost… Testosterone Testo… AR "" HMRbase… Non-prote…
2 Testosterone-Testost… Testosterone Testo… AR "" HMRbase… Non-prote…
3 Thyrostimulin-GPHB5_… TSH GPHB5… TSHR "" PMID: 1… Non-prote…
4 Thyrotropin-CGA_TSHB… TSH CGA_T… TSHR "" PMID: 1… Non-prote…
5 Triiodothyronine-T3-… Triiodothyr… T3-DI… THRA "" HMRbase… Non-prote…
6 Triiodothyronine-T3-… Triiodothyr… T3-DI… THRB "" HMRbase… Non-prote…Preprocess expression data
To infer the cell state-specific communications, CellChat identifies over-expressed ligands or receptors in one cell group and then identifies over-expressed ligand-receptor interactions if either ligand or receptor are over-expressed.
# subset the expression data of signaling genes for saving computation cost
cc_covid@DB <- CellChatDB
cc_covid <- subsetData(cc_covid) # This step is necessary even if using the whole database
future::plan("multisession", workers = 4) # do parallel
cc_covid <- identifyOverExpressedGenes(cc_covid)
cc_covid <- identifyOverExpressedInteractions(cc_covid)The number of highly variable ligand-receptor pairs used for signaling inference is 2205 # The number of highly variable ligand-receptor pairs used for signaling inference is 2212 Same for the flu
# subset the expression data of signaling genes for saving computation cost
cc_flu@DB <- CellChatDB
cc_flu <- subsetData(cc_flu) # This step is necessary even if using the whole database
future::plan("multisession", workers = 4) # do parallel
cc_flu <- identifyOverExpressedGenes(cc_flu)
cc_flu <- identifyOverExpressedInteractions(cc_flu)The number of highly variable ligand-receptor pairs used for signaling inference is 1720 # The number of highly variable ligand-receptor pairs used for signaling inference is 1728 Inference of cell-cell communication network
From the CellChat vignette: CellChat infers the biologically significant cell-cell communication by assigning each interaction with a probability value and peforming a permutation test. CellChat models the probability of cell-cell communication by integrating gene expression with prior known knowledge of the interactions between signaling ligands, receptors and their cofactors using the law of mass action.
The number of inferred ligand-receptor pairs clearly depends on the method for calculating the average gene expression per cell group. trimean approximates 25% truncated mean, implying that the average gene expression is zero if the percent of expressed cells in one group is less than 25%.
By setting this stringent parameters we ensure we are capturing robust signal and removing spurious cell-cell communication.
cc_covid <- computeCommunProb(cc_covid, type = "triMean")triMean is used for calculating the average gene expression per cell group.
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2024-03-08 10:32:05.470454]"
[1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2024-03-08 10:34:34.199726]"cc_flu <- computeCommunProb(cc_flu, type = "triMean")triMean is used for calculating the average gene expression per cell group.
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2024-03-08 10:34:36.527255]"
[1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2024-03-08 10:36:04.791565]"# If we wanted to be more lax we could use truncatedMean
# By setting trim - 0.1 we are only removing the top and bottom 10% expressing cells for that gene
# cc_covid <- computeCommunProb(cc_covid, type = "truncatedMean", trim = 0.1)
# cc_flu <- computeCommunProb(cc_flu, type = "truncatedMean", trim = 0.1)Summarize CCC to a signaling pathway level
We can summarize the individual ligand-receptor interactions into pathways for a broader and more functional look at the data.
cc_covid <- computeCommunProbPathway(cc_covid)
cc_flu <- computeCommunProbPathway(cc_flu)
## Look at all the pathways availavle
cc_covid@netP$pathways [1] "MHC-I" "MHC-II" "MIF" "CD99"
[5] "GALECTIN" "CypA" "CLEC" "ICAM"
[9] "ANNEXIN" "ADGRE" "CD45" "CCL"
[13] "APP" "Prostaglandin" "BAFF" "SELPLG"
[17] "LCK" "PECAM1" "CD23" "RESISTIN"
[21] "TGFb" "Cholesterol" "VISFATIN" "IL16"
[25] "SEMA4" "PARs" "GRN" "TNF"
[29] "CXCL" "GP1BA" "COLLAGEN" "CD40"
[33] "CysLTs" "CD160" "CD86" "LAIR1"
[37] "CD6" "PECAM2" "IFN-II" "ESAM"
[41] "ICOS" "BAG" "FLT3" "IL1"
[45] "TWEAK" "NECTIN" "MPZ" Calculate the aggregated cell-cell communication network by counting the number of links or summarizing the communication probability.
cc_covid <- aggregateNet(cc_covid)
cc_flu <- aggregateNet(cc_flu)Lastly we can compute the network centrality scores of the cell types which we can use to identify key senders, receivers, mediators…
# the slot 'netP' means the inferred intercellular communication network of signaling pathways
cc_covid <- netAnalysis_computeCentrality(cc_covid, slot.name = "netP")
cc_flu <- netAnalysis_computeCentrality(cc_flu, slot.name = "netP") Visualize inferred IFN-II signaling network
With Heatmap
object.list <- list(covid = cc_covid, flu = cc_flu)
par(mfrow = c(1, 2))
plt_ls <- lapply(seq_len(length(object.list)), function(i) {
netVisual_heatmap(
object.list[[i]],
signaling = "IFN-II",
color.heatmap = "Reds",
title.name = glue::glue("IFN-II signaling network - {names(object.list)[i]}"))
})
draw(plt_ls[[1]] + plt_ls[[2]], ht_gap = unit(.5, "cm"))
Merge Covid and Flue objects
Now we can merge the CellChat objects and work with them together. The rest of the analysis follows the comparison vignette.
cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)The cell barcodes in merged 'meta' is AAACCCAAGAATGTTG-12 AAACCCAAGGCCTGCT-12 AAACCCAAGGGCAATC-1 AAACCCAAGTAAGGGA-18 AAACCCACAAGAGCTG-17 AAACCCACACAGTCGC-20 Save objects in case we need them for later use
saveRDS(object.list, file = "../data/cellchat_ls.rds")
# object.list <- readRDS(file = "../data/cellchat_ls.rds")
# rm(object.list); gc()Differential interaction analysis
Let’s start by comparing the the total number of interactions and the sum of the weights
gg1 <- compareInteractions(
cellchat,
show.legend = FALSE,
group = c(1,2))
gg2 <- compareInteractions(
cellchat,
show.legend = FALSE,
group = c(1,2),
measure = "weight")
gg1 + gg2
With this plot we can visualize how the covid dataset has a slightly larger number of interactions but the interaction strength within the flu dataset are stronger/more likely.
We can also take a look at how these interactions are different between covid and flu. In the colorbar, red (or blue) represents increased (or decreased) signaling in the second dataset compared to the first one. In our case red is stronger in Flu and blue in covid
gg1 <- netVisual_heatmap(cellchat)
gg2 <- netVisual_heatmap(cellchat, measure = "weight")
gg1 + gg2
Next we can see how much the incoming and outgoing interaction signals are changing between conditions. This requires some data prep so we are showing the plot below:
slot.name <- "netP"
x.measure <- "outdeg"
y.measure <- "indeg"
signaling <- NULL
df_ls <- lapply(names(object.list), function(nm) {
object <- object.list[[nm]]
centr <- slot(object, slot.name)$centr
outgoing <- matrix(0, nrow = nlevels(object@idents), ncol = length(centr))
incoming <- matrix(0, nrow = nlevels(object@idents), ncol = length(centr))
dimnames(outgoing) <- list(levels(object@idents), names(centr))
dimnames(incoming) <- dimnames(outgoing)
for (i in 1:length(centr)) {
outgoing[, i] <- centr[[i]][[x.measure]]
incoming[, i] <- centr[[i]][[y.measure]]
}
outgoing.cells <- rowSums(outgoing)
incoming.cells <- rowSums(incoming)
num.link <- aggregateNet(object, signaling = signaling, return.object = FALSE,
remove.isolate = FALSE)$count
num.link <- rowSums(num.link) + colSums(num.link) - diag(num.link)
df <- data.frame(
outgoing.cells,
incoming.cells,
names(incoming.cells),
num.link,
nm)
colnames(df) <- c(paste0("outgoing_", nm), paste0("incoming_", nm), "labels", paste0("count_", nm))
df$labels <- factor(df$labels, levels = names(incoming.cells))
df
})Visualize how the interactions changed
require(gridExtra)
library(colorBlindness)
library(RColorBrewer)
# Define the number of colors you want
nb.cols <- length(unique(se$Celltype))
col_pal <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)
names(col_pal) <- sort(unique(se$Celltype))
# Create a ggplot with 18 colors
df_ls[[1]] %>%
left_join(df_ls[[2]], by = "labels") %>%
ggplot() +
geom_point(aes(x = outgoing_flu, y = incoming_flu, color = labels,
fill = "Flu", size = count_flu), alpha = 0.5) +
geom_point(aes(x = outgoing_covid, y = incoming_covid, color = labels,
fill = "Covid", size = count_covid), size = 3, pch = 21) +
geom_segment(aes(
x = outgoing_covid,
y = incoming_covid,
xend = outgoing_flu,
yend = incoming_flu,
color = labels),
arrow = arrow(angle = 40, length = grid::unit(.35, "cm")),
show.legend = FALSE) +
theme_classic() +
guides(
size = FALSE,
shape = FALSE,
colour = guide_legend(override.aes = list(size = 5))) +
scale_color_manual(values = col_pal) +
labs(
title = "Cell-Cell interaction strength shifts from Covid to Flu",
# x = "Outgoing strength",
# y = "Incoming strength",
y = NULL, x = NULL,
color = "Cell Type",
fill = "Disease") +
scale_fill_manual(values = c("white", "black")) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line()) +
scale_y_continuous(labels = function(x) {
case_when(x == 25 ~ "25 Incoming\nStrength",
TRUE ~ as.character(round(x)))
}) +
scale_x_continuous(labels = function(x) {
case_when(x == 10 ~ "10\nOutgoing\nStrength",
TRUE ~ as.character(round(x)))
}) +
coord_fixed()
In this case we can see how CD8 T cells greatly increase their incoming signal as well as intermediate monocytes. NKs and CD4 EM-like, in turn, greatly increase their outgoing signal while not increasing their incoming signal.
CellChat allows us to explore signalling changes at the pathway level within one specific population between both conditions. Positive values highlight increases in the second group (Flu) while negative values are increases in the first group (covid).
netAnalysis_signalingChanges_scatter(cellchat, idents.use = "CD4, EM-like")
Differential pathway signalling
Look at broad changes in pathway signalling across the entire dataset
gg1 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = F, do.stat = TRUE)
gg1 + gg2
Diving into a pathway
We also need to take a look at which are the LR interactions driving a pathway. We can easily look at this using the function XX
(netAnalysis_contribution(object.list[[1]], signaling = "MHC-II") +
labs(title = "Contributions of each L-R pair in Covid")) /
(netAnalysis_contribution(object.list[[2]], signaling = "MHC-II") +
labs(title = "Contributions of each L-R pair in Flu"))
Comparing all the pathways at once
Lastly, we can also take a look at all the pathway changes at once by cell type between both conditions.
library(ComplexHeatmap)
# combining all the identified signaling pathways from different datasets
pathway.union <- union(object.list[["covid"]]@netP$pathways, object.list[["flu"]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(
object.list[["covid"]],
pattern = "outgoing",
signaling = pathway.union,
title = "covid",
width = 5,
height = 12)
ht2 <- netAnalysis_signalingRole_heatmap(
object.list[["flu"]],
pattern = "outgoing",
signaling = pathway.union,
title = "flu",
width = 5,
height = 12)
draw(ht1 + ht2, ht_gap = unit(.5, "cm"))
# combining all the identified signaling pathways from different datasets
pathway.union <- union(object.list[["covid"]]@netP$pathways, object.list[["flu"]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(
object.list[["covid"]],
pattern = "incoming",
signaling = pathway.union,
title = "covid",
width = 5,
height = 12)
ht2 <- netAnalysis_signalingRole_heatmap(
object.list[["flu"]],
pattern = "incoming",
signaling = pathway.union,
title = "flu",
width = 5,
height = 12)
draw(ht1 + ht2, ht_gap = unit(.5, "cm"))
Visualizing L-R interactions between cell types
levels(cellchat@idents$joint) [1] "B cell, IgG-" "B cell, IgG+" "CD4, EM-like"
[4] "CD4, non-EM-like" "CD8, EM-like" "CD8, non-EM-like"
[7] "classical Monocyte" "DC" "intermediate Monocyte"
[10] "NK cell" "nonclassical Monocyte" "Platelet"
[13] "RBC" netVisual_bubble(
cellchat,
sources.use = "DC",
targets.use = c("NK cell", "CD4, EM-like"),
comparison = c(1, 2),
angle.x = 45)
We can see here how the MHC-II signalling is present between DCs & CD4 and completely absent with NKs.
Session Info
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] RColorBrewer_1.1-3 colorBlindness_0.1.9 gridExtra_2.3
[4] ComplexHeatmap_2.16.0 CellChat_2.1.2 Biobase_2.60.0
[7] BiocGenerics_0.46.0 ggplot2_3.5.0 igraph_2.0.2
[10] dplyr_1.1.4 Seurat_5.0.2 SeuratObject_5.0.1
[13] sp_2.1-3
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.22 splines_4.3.1 later_1.3.2
[4] tibble_3.2.1 polyclip_1.10-6 ggnetwork_0.5.13
[7] fastDummies_1.7.3 lifecycle_1.0.4 rstatix_0.7.2
[10] doParallel_1.0.17 vroom_1.6.3 globals_0.16.2
[13] processx_3.8.2 lattice_0.21-8 MASS_7.3-60
[16] backports_1.4.1 magrittr_2.0.3 sass_0.4.8
[19] plotly_4.10.4 rmarkdown_2.25 jquerylib_0.1.4
[22] yaml_2.3.8 remotes_2.4.2.1 httpuv_1.6.14
[25] NMF_0.26 sctransform_0.4.1 spam_2.10-0
[28] sessioninfo_1.2.2 pkgbuild_1.4.2 spatstat.sparse_3.0-3
[31] reticulate_1.35.0.9000 cowplot_1.1.3 pbapply_1.7-2
[34] abind_1.4-5 pkgload_1.3.2.1 Rtsne_0.17
[37] presto_1.0.0 purrr_1.0.2 circlize_0.4.15
[40] IRanges_2.34.1 S4Vectors_0.38.1 ggrepel_0.9.5
[43] irlba_2.3.5.1 listenv_0.9.1 spatstat.utils_3.0-4
[46] goftest_1.2-3 RSpectra_0.16-1 spatstat.random_3.2-2
[49] fitdistrplus_1.1-11 parallelly_1.37.0 svglite_2.1.3
[52] leiden_0.4.3.1 codetools_0.2-19 tidyselect_1.2.0
[55] shape_1.4.6 farver_2.1.1 matrixStats_1.2.0
[58] stats4_4.3.1 spatstat.explore_3.2-6 jsonlite_1.8.8
[61] BiocNeighbors_1.18.0 GetoptLong_1.0.5 ellipsis_0.3.2
[64] progressr_0.14.0 ggalluvial_0.12.5 ggridges_0.5.6
[67] survival_3.5-7 iterators_1.0.14 systemfonts_1.0.4
[70] foreach_1.5.2 tools_4.3.1 sna_2.7-2
[73] ica_1.0-3 Rcpp_1.0.12 glue_1.7.0
[76] xfun_0.42 usethis_2.2.2 withr_3.0.0
[79] BiocManager_1.30.22 fastmap_1.1.1 fansi_1.0.6
[82] callr_3.7.3 digest_0.6.34 gridGraphics_0.5-1
[85] R6_2.5.1 mime_0.12 colorspace_2.1-0
[88] Cairo_1.6-1 scattermore_1.2 tensor_1.5
[91] spatstat.data_3.0-4 utf8_1.2.4 tidyr_1.3.1
[94] generics_0.1.3 data.table_1.15.0 FNN_1.1.4
[97] prettyunits_1.1.1 httr_1.4.7 htmlwidgets_1.6.4
[100] uwot_0.1.16 pkgconfig_2.0.3 gtable_0.3.4
[103] registry_0.5-1 lmtest_0.9-40 htmltools_0.5.7
[106] carData_3.0-5 profvis_0.3.8 dotCall64_1.1-1
[109] clue_0.3-64 tidyverse_2.0.0 scales_1.3.0
[112] png_0.1-8 knitr_1.45 rstudioapi_0.15.0
[115] tzdb_0.4.0 reshape2_1.4.4 rjson_0.2.21
[118] coda_0.19-4.1 statnet.common_4.9.0 nlme_3.1-163
[121] curl_5.2.0 cachem_1.0.8 zoo_1.8-12
[124] GlobalOptions_0.1.2 stringr_1.5.1 KernSmooth_2.23-22
[127] parallel_4.3.1 miniUI_0.1.1.1 pillar_1.9.0
[130] vctrs_0.6.5 RANN_2.6.1 ggpubr_0.6.0
[133] urlchecker_1.0.1 promises_1.2.1 car_3.1-2
[136] xtable_1.8-4 cluster_2.1.4 evaluate_0.23
[139] magick_2.7.5 readr_2.1.4 cli_3.6.2
[142] compiler_4.3.1 rlang_1.1.3 crayon_1.5.2
[145] rngtools_1.5.2 ggsignif_0.6.4 future.apply_1.11.1
[148] labeling_0.4.3 ps_1.7.5 plyr_1.8.9
[151] fs_1.6.3 stringi_1.8.3 network_1.18.2
[154] BiocParallel_1.34.2 viridisLite_0.4.2 deldir_2.0-2
[157] gridBase_0.4-7 munsell_0.5.0 lazyeval_0.2.2
[160] devtools_2.4.5 spatstat.geom_3.2-8 Matrix_1.6-5
[163] RcppHNSW_0.6.0 hms_1.1.3 patchwork_1.2.0
[166] bit64_4.0.5 future_1.33.1 shiny_1.8.0
[169] ROCR_1.0-11 broom_1.0.5 memoise_2.0.1
[172] bslib_0.6.1 bit_4.0.5